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We use a theoretical frame-work to analytically assess temporal prediction error functions on von-Karman 
turbulence when a zonal representation of wave-fronts is assumed. Linear prediction models analysed include 
auto-regressive of order up to three, bilinear interpolation functions and a minimum mean square error predic¬ 
tor.This is an extension of the authors’ previously published work Q in which the efficacy of various temporal 
prediction models was established. Here we examine the tolerance of these algorithms to specific forms of 
model errors, thus defining the expected change in behaviour of the previous results under less ideal condi¬ 
tions. Results show that ±100% wind-speed error and ±50 deg are tolerable before the best linear predictor 
delivers poorer performance than the no-prediction case. 
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Temporal prediction of the atmosphere is a much de¬ 
bated topic. The purpose of prediction is to reduce the 
error due to servo lag since the turbulence profile changes 
rapidly (on timescales of a few milliseconds) during the 
time that it takes to gather sensor information and to 
compute corrections. 

Unlike single-conjugated AO systems, in tomographic 
AO direct access to an estimate of the layered wave- 
front (WF) is provided at the end of the tomographic 
step. With turbulence estimated in a discrete number of 
layers, the frozen-flow approximation can now be called 
upon with a higher degree of fidelity [l| . 

Temporal prediction is useful for both Multi- 
Conjugate AO (MCAO) and Multi-Object AO (MOAO) 
as a means to increase the sky-coverage ii- Allowing 
for greater integration times whilst compensating for the 
lag error by applying a predictive algorithm enables the 
system to guide on fainter sources [2|. Provided infor¬ 
mation on the dynamics of the atmospheric turbulence 
is available or can be construed, one should be able to 
obtain a more accurate estimate of the WF at the time 
a set of commands is applied to the deformable mirror 
(DM) and therefore improve performance. Lag errors 
are also considered a serious limitation in high-contrast 
imaging systems where the broadening of the PSF along 
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the main axis of wind-blown turbulence severely limits 
contrast at small separations [^. 

In this letter we discuss alternatives for time- 
progressing the atmospheric wave-fronts namely a near- 
Markovian model, auto-regressive models and spatial 
shifting under frozen-flow E S The main goal of 
our work is to provide insight into the accuracy and ro¬ 
bustness bounds of such models in view of a contem¬ 
porary application to the Raven science and technology 
demonstrator installed on the Subaru Telescope Q. 

Under the hypothesis that the turbulent atmosphere 
is a sum of L thin layers located at a discrete number of 
different altitudes hi, the aperture-plane phase 4>{p, 0, t) 
indexed by the bi-dimensional spatial coordinate vector 
p = [px, Py\ in direction 6 = [9x, 9y\ at time t is defined 
as 

L 

4>ip,e,t) = U0(p{p,t) = '^u}ppi{p-ir hie,t) (1) 

1=1 

where ipi{p,t) is the /‘^-layer wave-front, U!i is the 
layer strength and Hg is a propagation operator in the 
near-held approximation that relates the aperture-plane 
WF to the wave-front dehned over a discrete number 
of L layers in the volume by adding and interpolating 
tpi{p-\-hi9,t) on the aperture-plane computational grid. 

In the following we assume a point-wise representa¬ 
tion of the WF. In both MCAO and MOAO, pseudo 
open-loop or straight open-loop estimation of the phase 
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is performed from noisy measurements involving derived 
quantities related to through the operator G 

with additive noise r]{p,t) in Eq. ([2|) - by means of 
a linear wave-front reconstructor. The goal is thus to 
estimate 


^ip,t + T) =W[Gcj){p,t) + T]{p,t)] ( 2 ) 

where W = WtWp is a linear reconstructor that per¬ 
forms spatial (indicated by subscript p) and temporal 
(indicated by subscript t) phase estimation sequentially. 

The progression of WFs with time is not a well-defined 
process. We use the simplifying Taylor or frozen-flow 
approximation; under this hypothesis, time-progressing 
the phase is equivalent to spatially shifting it in both 
’x’ and ’y’ directions by [Aa;,Ay] = [—VxT^—VyT], the 
product of the components of the wind velocity vector, 
V = [vx'jVy] with the lag, r. Thus the turbulence profile 
at time (t + r) and height I is obtained simply by trans¬ 
lating the profile by Ap = vr. The temporally shifted 
phase can then be expressed in terms of the spatial shift, 
ip{p,t + T) = ip{p — \'T,t)- Under this assumption, linear 
temporal predictive models, A{t), are discussed next. 

Suppose the current and future WFs are related by 

ip{p,t + T) = A{T)ip{p,t) (3) 

where A(t) is a block-diagonal linear operator (atmo¬ 
spheric layers are considered independent) which time- 
progresses the WF. For the sake of compactness the layer 
dependence is dropped in the following but is subtended. 

The linear model that minimizes the mean square pre¬ 
diction residual phase variance is @ 

A* = aTgimn (^\\^p{t + t) - A*, ( 4 ) 


where L 2 stands for Euclidean norm, is the telescope 
aperture and (•) means ensemble averaging over the tur¬ 
bulence statistics. It is best to introduce the notation 
ip{t + t) A and ip{t) A which are the sampled 

phase vectors with r lag in between. 

The solution to Eq. O is readily found to be the 
minimum mean square error (MMSE) estimator 

= {‘Pk+i<fik) ^, (5) 

where i{p) = VkiiP ~ v;t). This predictor is, in 
what follows, called the Spatio-Angular predietor (SA). 
The covariance matrices are computed by sampling the 
covariance function with distances corresponding to all 
baselines between phase points in the aperture 


c^{p) 


r(ii/6) 

\ro) 25/6^8/3 




( 6 ) 


where Lq is the outer scale of turbulence, rp Fried’s 
parameter, F the ’gamma’ function and finally ATs/p a 


modified Bessel function of the third order. Parameter 
p = IpI is the distance between any two phase points in 
the bi-dimensional plane. 

The operator A* is invertible. With the turbu¬ 
lence covariance matrix 51,^ translation invariant, then 

Furthermore, Eq. is a general method to generate 
and predict phase in a 2D plane (any wind velocity can 
be used) according to the Markovian model 

^k+i=A*^k + ^l ( 7 ) 

where is an excitation noise whose covariance func¬ 
tion S* is fixed to guarantee proper turbulence statis¬ 
tics from either a Kolmogorov or von-Karman model. 
Hence S* = - A*i:^A*’'^, since {<Pk+ii‘Pk+iV) = 

+ s:._Note A 

which, due to stationarity, loses it temporal dependence. 

However, it has previously been pointed out @ that 
the phase over the aperture is not Markovian since some 
information on the wind-blown portion that leaves the 
telescope aperture is dropped. We will nevertheless use 
this near-Markovian model since the temporal update 
rate of the phase estimates is often of the order of, or 
below, the coherence time of the turbulence. 

The near-Markovian model of Eq. 0 generalizes 
the one from Q for simulating infinitely long, non¬ 
stationary phase screens to fractional pixel shifts in a 
bi-dimensional plane. In Assemat’s method, a sub-set 
of the columns of A* is used (most commonly 2) to in¬ 
clude bounded-region correlations. 

At this point, two important observations motivate 
the discussion that follows: i) A* is strongly diagonally 
dominated, suggesting that simpler diagonal models, i.e. 
point-wise, could potentially be applied - and indeed 
they have been extensively used in AO simulations in 
the form of auto-regressive (AR) models @ albeit with 
Zernike functions; ii) A* works out to shift the phase 
screen in the appropriate wind direction with shifts given 
by |v/|cos(0i) and \vi\sin{9i) for layer I = 1 • ■ • L. Fur¬ 
thermore, the new turbulence that enters in the tele¬ 
scope aperture is estimated using spatial correlations 
with all the points in the aperture - see Fig. [21 

As mentioned above, this model is specific for frozen- 
flow, which is but a part of the actual atmospheric dis¬ 
turbances [I| . The robustness to parameter mismatches 
of models with such an underlying assumption must 
therefore be well understood. 

An AR model of order n is defined by the recursion 

‘Pk+I = f {^kr ■ ■ Vk-n-i) A ( 8 ) 

where /(■••) is a linear function yet to be defined and 
£k is a Gaussian-distributed, spectrally white, zero-mean 
random sequence with covariance given in ref. 0. 

Diagonal AR models are suitable for systems con¬ 
strained by real-time computational burden. Using sim¬ 
pler diagonal auto-regressive models can be quite ap¬ 
pealing as they circumvent A* being a dense matrix. 
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Although these relatively coarse models are not adapted 
to simulating atmospheric turbulence, they are used in¬ 
stead for prediction when embedded in the reconstruc¬ 
tor and plugged into the dynamical controller, as is done 
with Kalman filtering, for off-line computation of opti¬ 
mal gains. The AR models are included here because 
they have been shown to provide stable and improved 
error reduction within that framework ai- 

Consider the following models 

Vk+i = ^ARiVk + er' (9) 

<Pk+l = -^AR2<^k + ^AR2<^k-l + ^k^^ 

<Pk+l = ^ARsVk + BARaV’fc-l + CAR3¥’fe-2 + 

Figure [T] depicts the temporal auto-correlation func¬ 
tions for the AR models with the fitted second order 
model achieving an overall best fit. 

A practical alternative to model fitting is the use of the 
Yule-Walker equations for auto-regressive model compu¬ 
tation to identify ARn models from a discretized covari¬ 
ance function at the AO system’s sampling frequency - 
also shown in Fig. [T] 


Temporal Auto-Correlation Functions 



Fig. 1. Left: Temporal auto-correlation fmrctions for tur¬ 
bulence with To — 0.155 m, Lq = 25 m and identified ARl, 
AR2 and AR3 models. The vertical dotted line indicates the 
fittiirg section. 

Figure [2] depicts an analytical covariance matrix and a 
thresholded Spatio-Angular prediction matirx. A* from 
Eq. ([5]). It is apparent that the optimal one-step esti¬ 
mator is largely dominated by a diagonal term that sim¬ 
ply shifts the phase in the appropriate direction. Note 
however that points outside the pupil entering the tele¬ 
scope (upper band in Figl^right) are computed using 
their respective covariance function with points inside 
the pupil. We thus expect the bilinear interpolation to 
be a reasonable approximation to the Spatio-Angular 
model but with degraded performance. It will be tested 
and assessed in numerical simulations. 

Consider the general formulation for the temporal lag 
error that is compatible with the case of atmospheric 



Fig. 2. Left: Spatial phase covariance matrix (single-layer, 
squared pupil, sampled with 33 points across). Right: Matrix 
A* for an horizontal shift of 0.5m. Sampling 0.25m. A* is 
largely dominated by a diagonal term that simply shifts the 
phase in the appropriate direction. The same stands when 
the shifts are not integer factors of the spatial sampling. In 
that case. A* becomes more populated with up to four main 
diagonals that interpolate every single point based on the 
neighbouring 4 points. Matrix shown for entries > 10“"^. 

prediction 

= (||He {<fik+i - Vfc+i) 111,(0)) (10) 

where is a layered phase estimate at time step k. 

In the no prediction case, the estimated phase is sim¬ 
ply a replication of the phase at the previous time step, 
iPk = ‘fik-i- The temporal lag error from Ea. dTUl) be¬ 
comes 

al{T) = trace {He HJ} 4 (n) 

which is simply the temporal structure function of phase, 
with the 1-step covariance matrix s; = s^(p' = p + 
|v|r) computed from Eq. ®. 

These temporal structure functions can now be ex¬ 
panded for the case of predicted phase. Developing Eq. 
(fTUl) using the measurement models given in Eq. [5]yields, 

1. First order models (ARl or SA) 

crl^i{r) = trace (He (S^ + HJ} , 

( 12 ) 

2. Second order models (AR2) 

f^p= 2 ('^) = trace {Hq {'^ip + AS,^A -f 

-2AS; + 2AS;bT - 2BS") HJ} (13) 

with S" = S<p(p" = p + 2|v|t). 

3. Third order models (AR3) 

<^p =3 (t) = trace {He (Sy, -f- AS,^A -f- BS,^B-t- 
CS^C - 2AS; - 2BS" - 2CS"' 
+2AS;bT + 2BS;c'^ + 2AS"CT) Hi} 

(14) 

with S"' = Sip(/o'" = p + 3|v|r). 
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These equations are revised from those given in 0, 
which contain a minor transcription error; the numerical 
results were not affected. The AR models perform very 
poorly when prediction is considered; in Fig. [3] the the¬ 
oretical temporal structure functions indicate that little 
improvement is gained from the no-prediction case when 
the AR models are considered. The SA model performs 
by far best as previously noted in Qj. 


Temporal Structure Functions 



Fig. 3. Theoretical temporal structure functions quoted in 
rad^@0.5fim under same conditions as Fig. [T] 

With the expectation that it is a reasonable approx¬ 
imation to the SA model, we built the bilinear spline 
interpolation model with two versions. In the first, new 
points entering the aperture are interpolated with no fur¬ 
ther information from within the pupil, thus assuming 
zeroed wave-front outside the pupil. In the second, the 
new rows are estimated from a correlation function, i. e. 
those new rows in A (Fig. [2|) are kept, the remainder 
being the two-dimensional bilinear spline weights. The 
results in Fig. |l](top) show that if one does not estimate 
the new points from a correlation function, the loss in 
performance is quite drastic, with the bilinear spline in¬ 
terpolator achieving best performance for a wind-speed 
model that is 1 /2 the true wind-speed. When the values 
of the new points are computed keeping the correspond¬ 
ing columns of A* then the performance enhancement is 
dramatic with a difference of some tens of nm rms. 

Also shown in Fig |T] are the analytical temporal error 
functions for a model mismatch in terms of wind-speed 
and wind direction for the ARl model and for the full 
MMSE predictor. Results for the ARl model show that, 
although it will never reduce temporal lag error, it does 
not increase it either even in very poorly identified con¬ 
ditions. Results for the MMSE model (and similarly 
for the bilinear interpolation with optimized new entries 
into the pupil) indicate that the wind speed of a given 
layer can be over estimated by up to a factor of 2 before 
lag error reaches the equivalent of no prediction. An er¬ 
ror in wind direction of up to 50° still reduces temporal 
lag error; results are given for a 50Hz sample rate, but 
the bounds hold in the range of lOHz - IkHz. 


Predictive Models Robustness, d=0.5m, x=0.02s, |Vg|=10m/s 



-0.5 0 0.5 1 1.5 2 

mismatch factor k, ]v|=k’|VQ| 


Predictive Models Robustness, d=0.5m, x=0.02s, |Vg|=10m/s 



Fig. 4. (Top:) Robustness wrt wind-speed. (Bottom:) Ro¬ 
bustness wrt wind direction. 
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